clear all


cd "$revdta"
set more off
 

***********PARTICIPATION MOMENTS BEFORE AND AFTER AGE 45 ****************************************************************************************************************
foreach type in $type {
local typename $typename
use ".\Full_Panel_2b.dta", clear
keep if sample_analytic==1 & educ <= 2

*drop people who are ever on SSI and not on DI
*bysort newid: egen everSSI=max(tr_ssi >0 & tr_ssi != . & DI ==0)
*bysort newid: egen everboth=max(tr_ssi >0 & DI==1)

*drop everSSI

replace DI = anyDI

gen anywork = hours > 0 & hours != .
gen emp_part = hours > 0 & hours <= 1500 & hours != . 
gen emp = hours > 1500 & hours != . 

*bysort newid: egen everworkhealthy=max(emp == 1 & DS == 0)
*keep if everworkhealthy == 1 & everSSI == 0
*drop everworkhealthy

bysort newid: egen firstage = min(age)
bysort newid: egen firsthealth = max(DS * (age==firstage) - 999*(age!=firstage))

*local var firsthealth
*local var chd_maxhealth
local var
if("`var'"== "firsthealth") {
 keep if firsthealth ==0
 tokenize `""fh""'
 }
if("`var'" == "chd_maxhealth") {
keep if chd_maxhealth>=1 & chd_maxhealth<=3 
tokenize `""chd""'
}
*I lose most people before 1996 from child health, another 30% from restricting 
*to in-sample heads
*total: 2k households


gen sp_dis = sp_DS==1 | sp_DS==2

gen part_full = lw!=. & hours >=1500 
gen part_part = lw !=. & hours < 1500 & hours >=500

gen sp_part_full = sp_lw !=. & sp_hours >=1500
gen sp_part_part = sp_lw != . & sp_hours < 1500 & sp_hours >=500

egen everwork = max(part_full), by(newid)


egen part_int = group(part_full sp_part_full)
egen part_intpart = group(part_part sp_part_part)
replace sp_DS = sp_DS >=1
gen     age_band=1 if age<=45
replace age_band=2 if age>45

gen     sp_age_band=1 if agew<=45
replace sp_age_band=2 if agew>45 & agew != .
tab age_band,gen(ad)


merge m:1 newid using "./intermediary/sp_temp_types"
keep if _merge ==3
ren type_fix type_lp
preserve
ren `type' type_fix


/*
replace food=. if food>1199987		
replace fout=. if food>129999
egen totfood=rsum(food fout)
replace rent=0.06*house if rent==0 & house>0
egen spending=rsum(totfood health educexp util childcare transport rent homeins)
replace spending=spending/price
*/
replace spending = spending/(1+0.5*(fsize-kids-1)+0.3*kids)
replace spending = . if year <= 1996

keep spending age newid type_fix
gen obs = 1
tempfile times
save `times', replace
drop spending
reshape wide obs, i(newid) j(age) 
drop newid
order type_fix, first
foreach x of varlist obs* {
replace `x' = 0 if `x' == .
}
forvalues j=63(1)72{
gen obs`j'=0
}
save "$out\psid_obs`1'`typename'.dta", replace

use `times', clear
keep if spending != .
drop spending
reshape wide obs, i(newid) j(age)
drop newid
order type_fix, first
foreach x of varlist obs* {
replace `x' = 0 if `x' == .
}
forvalues j=63(1)72{
gen obs`j'=0
}
save "$out\psid_obs_cons`1'`typename'.dta", replace

restore
}
ren $type type_fix

preserve

gen sp_birthybin = round(sp_birthy,5)
gen n = 1
collapse (mean) sp_part_full (count) n, by(agew sp_birthybin)
keep if n >=20
levelsof sp_birthybin, local(bins)
local plot 
local p = 1
foreach x of local bins {
local plot `plot' (connected sp_part agew if sp_birthybin==`x', legend(label(`p' "Bin `x'")))
local p = `p' +1
}

twoway `plot'
restore
***** Average participation by age ****
matrix b=0
matrix V=0

cap program drop doit
program def doit
foreach suffix in full part {
         local i=0
	   while `i' < 3 {
		   local j=1	
		   while `j' < 3 {
				reg part_`suffix' if DS==`i' & age_band==`j' ,vce(cluster newid)
				matrix b=b\e(b)
				matrix V=V\e(V)
	         local j=`j'+1
                 }
         local i=`i'+1
                 }
				 }
				 
end

cap program drop sp_doit
program def sp_doit

foreach suffix in full part {
         local i=0
	   while `i' < 2 {
		   local j=1	
		   while `j' < 3 {
				reg sp_part_`suffix' if sp_DS==`i' & sp_age_band==`j'  ,vce(cluster newid)
				matrix b=b\e(b)
				matrix V=V\e(V)
	         local j=`j'+1
                 }
         local i=`i'+1
                 }
				 }
end

qui doit
matrix VV=0
forvalues j=2(1)13	{
	matrix temp=sqrt(V[`j',1])
	matrix VV=VV\temp
	}
matrix avgpart_def1=b[2..13,1],(VV[2..13,1])

matrix rownames avgpart_def1 = Pr(E=2/L=0,Age<=45) Pr(E=2/L=0,Age>45) 	///
								 Pr(E=2/L=1,Age<=45) Pr(E=2/L=1,Age>45) ///
								 Pr(E=2/L=2,Age<=45) Pr(E=2/L=2,Age>45) ///
								 Pr(E=1/L=0,Age<=45) Pr(E=1/L=0,Age>45) 	///
								 Pr(E=1/L=1,Age<=45) Pr(E=1/L=1,Age>45) ///
								 Pr(E=1/L=2,Age<=45) Pr(E=1/L=2,Age>45) 
matrix colnames avgpart_def1 = EST SE	
***Table 4 Panel B
di "marrital: all"
matlist avgpart_def1, title(Empl.rate(Year t),DS(Intvw. time t)) cspec(o4& %25s | %5.4f & %5.4f & ) rspec(&-&&&&&&&&&&&&)
mat2txt, matrix(avgpart_def1) saving("$out/workmoments`1'_agg") replace


foreach m in 0 1 {
preserve
keep if married == `m'
qui doit
matrix VV=0
forvalues j=2(1)13	{
	matrix temp=sqrt(V[`j',1])
	matrix VV=VV\temp
	}
matrix avgpart_def1=b[2..13,1],(VV[2..13,1])

matrix rownames avgpart_def1 = Pr(E=2/L=0,Age<=45) Pr(E=2/L=0,Age>45) 	///
								 Pr(E=2/L=1,Age<=45) Pr(E=2/L=1,Age>45) ///
								 Pr(E=2/L=2,Age<=45) Pr(E=2/L=2,Age>45) ///
								 Pr(E=1/L=0,Age<=45) Pr(E=1/L=0,Age>45) 	///
								 Pr(E=1/L=1,Age<=45) Pr(E=1/L=1,Age>45) ///
								 Pr(E=1/L=2,Age<=45) Pr(E=1/L=2,Age>45) 
matrix colnames avgpart_def1 = EST SE	
***Table 4 Panel B
di "marrital: `m'"
matlist avgpart_def1, title(Empl.rate(Year t),DS(Intvw. time t)) cspec(o4& %25s | %5.4f & %5.4f & ) rspec(&-&&&&&&&&&&&&)
mat2txt, matrix(avgpart_def1) saving("$out/workmoments_`m'`1'") replace

matrix b = 0
matrix V = 0
qui sp_doit

matrix VV=0
forvalues j=2(1)9	{
	matrix temp=sqrt(V[`j',1])
	matrix VV=VV\temp
	}
matrix avgpart_def1=b[2..9,1],(VV[2..9,1])

matrix rownames avgpart_def1 = Pr(E=2/L=0,Age<=45) Pr(E=2/L=0,Age>45) 	///
								 Pr(E=2/L=1,Age<=45) Pr(E=2/L=1,Age>45) ///
								 Pr(E=1/L=0,Age<=45) Pr(E=1/L=0,Age>45) 	///
								 Pr(E=1/L=1,Age<=45) Pr(E=1/L=1,Age>45) 
matrix colnames avgpart_def1 = EST SE	
***SPOUSAL: Table 4 Panel B
di "marrital: `m'"
matlist avgpart_def1, title(Empl.rate(Year t),DS(Intvw. time t)) cspec(o4& %25s | %5.4f & %5.4f & ) rspec(&-&&&&&&&&)
mat2txt, matrix(avgpart_def1) saving("$out/spworkmoments_`m'`1'") replace

matrix b = 0
matrix V = 0

restore
}

*joint distribution of work, at working ages:
***** Average participation by age ****
matrix b=0
matrix V=0

cap program drop doit
program def doit
foreach suffix in 1 2 3 4  {
				reg jointpart`suffix' if age < 62 ,vce(cluster newid)
				matrix b=b\e(b)
				matrix V=V\e(V)
				 }
end

tab part_int, gen(jointpart)


matrix b = 0
matrix V = 0
qui doit
matrix VV=0
forvalues j=2(1)5	{
	matrix temp=sqrt(V[`j',1])
	matrix VV=VV\temp
	}
matrix avgpart_def1=b[2..5,1],(VV[2..5,1])

matrix rownames avgpart_def1 = Neither_Work Spouse_Works 	///
								 Head_Works Both_Work 
								 
matrix colnames avgpart_def1 = EST SE	
***Table 4 Panel B
matlist avgpart_def1, title(Empl.rate(Year t),DS(Intvw. time t)) cspec(o4& %25s | %5.4f & %5.4f & ) rspec(&-&&&&)
mat2txt, matrix(avgpart_def1) saving("$out/jointworkmoments`1'") replace

drop jointpart*

tab part_intpart, gen(jointpart)

matrix b = 0
matrix V = 0
qui doit
matrix VV=0
forvalues j=2(1)5	{
	matrix temp=sqrt(V[`j',1])
	matrix VV=VV\temp
	}
matrix avgpart_def1=b[2..5,1],(VV[2..5,1])

matrix rownames avgpart_def1 = Neither_Work Spouse_Works 	///
								 Head_Works Both_Work 
								 
matrix colnames avgpart_def1 = EST SE	
***Table 4 Panel B
matlist avgpart_def1, title(Empl.rate(Year t),DS(Intvw. time t)) cspec(o4& %25s | %5.4f & %5.4f & ) rspec(&-&&&&)
mat2txt, matrix(avgpart_def1) saving("$out/jointworkpartmoments`1'") replace


tab part_int if  age <= 62

*bysort sp_DS DS age_band: sum part sp_part 


********Stocks on DI by health, marital status, and age
gen old=age>45
sort newid age
bysort newid: gen worstDS = DS[1]
bysort newid: replace worstDS = max(DS[_n],worstDS[_n-1]) if _n > 1

matrix b=0
matrix V=0
cap program drop doit
program def doit
         local i=0
	   while `i' < 2 {
		   local j=0	
		   while `j' < 3 {
				reg DI if DS==`j' & old==`i', vce(cluster newid)
				matrix b=b\e(b)
				matrix V=V\e(V)
	         local j=`j'+1
                 }
         local i=`i'+1
                 }
end
qui doit
matrix VV=0
forvalues j=2(1)7	{
	matrix temp=sqrt(V[`j',1])
	matrix VV=VV\temp
	}
matrix stockDI_def1=b[2..7,1],VV[2..7,1]


matrix rownames stockDI_def1 =	Pr(DI=1/L=0,Age<=45) Pr(DI=1/L=1,Age<=45) Pr(DI=1/L=2,Age<=45) ///
									Pr(DI=1/L=0,Age>45) Pr(DI=1/L=1,Age>45) Pr(DI=1/L=2,Age>45)
matrix colnames stockDI_def1 = EST SE	

***Table 4 Panel C
matlist stockDI_def1,cspec(o4& %25s | %5.4f & %5.4f & ) rspec(&-&&&&&&) title(DI rate(Year t),DS(Intvw. time t+1))
mat2txt, matrix(stockDI_def1) saving("$out/stockDImoments`1'_agg") replace



matrix b=0
matrix V=0
cap program drop doit
program def doit
		 local m = 0
		 while `m' < 2 {
         local i=0
	   while `i' < 2 {
		   local j=0	
		   while `j' < 3 {
				reg DI if DS==`j' & old==`i' & married==`m', vce(cluster newid)
				matrix b=b\e(b)
				matrix V=V\e(V)
	         local j=`j'+1
                 }
         local i=`i'+1
                 }
				 local m = `m' + 1
				 }
				 local l = `l' + 1
end
qui doit
matrix VV=0
forvalues j=2(1)13	{
	matrix temp=sqrt(V[`j',1])
	matrix VV=VV\temp
	}
matrix stockDI_def1=b[2..13,1],VV[2..13,1]


matrix rownames stockDI_def1 =	Pr(DI=1/L=0,Age<=45,single) Pr(DI=1/L=1,Age<=45,single) Pr(DI=1/L=2,Age<=45,single) ///
									Pr(DI=1/L=0,Age>45,single) Pr(DI=1/L=1,Age>45,single) Pr(DI=1/L=2,Age>45,single) ///
									Pr(DI=1/L=0,Age<=45,married) Pr(DI=1/L=1,Age<=45,married) Pr(DI=1/L=2,Age<=45,married) ///
									Pr(DI=1/L=0,Age>45,married) Pr(DI=1/L=1,Age>45,married) Pr(DI=1/L=2,Age>45,married)
matrix colnames stockDI_def1 = EST SE	

***Table 4 Panel C
matlist stockDI_def1,cspec(o4& %25s | %5.4f & %5.4f & ) rspec(&-&&&&&&&&&&&&) title(DI rate(Year t),DS(Intvw. time t+1))
mat2txt, matrix(stockDI_def1) saving("$out/stockDImoments`1'") replace


********Stocks on DI by employment, health, marital status, and age

gen DI1 = DI == 1 & empnow==1
gen DI0 = DI == 1 & empnow==0

matrix b=0
matrix V=0
cap program drop doit
program def doit
	local l = 0
	while `l' < 2 {
		 local m = 0
		 while `m' < 2 {
         local i=0
	   while `i' < 2 {
		   local j=0	
		   while `j' < 3 {
				reg DI`l' if DS==`j' & old==`i' & married==`m' , vce(cluster newid)
				matrix b=b\e(b)
				matrix V=V\e(V)
	         local j=`j'+1
                 }
         local i=`i'+1
                 }
				 local m = `m' + 1
				 }
				 local l = `l'+1
		}
end
qui doit
matrix VV=0
forvalues j=2(1)25	{
	matrix temp=sqrt(V[`j',1])
	matrix VV=VV\temp
	}
matrix stockDI_def1=b[2..25,1],VV[2..25,1]


matrix rownames stockDI_def1 =	Pr(DI=1&W=0/L=0,Age<=45,single) Pr(DI=1&W=0/L=1,Age<=45,single) Pr(DI=1&W=0/L=2,Age<=45,single) ///
									Pr(DI=1&W=0/L=0,Age>45,single) Pr(DI=1&W=0/L=1,Age>45,single) Pr(DI=1&W=0/L=2,Age>45,single) ///
									Pr(DI=1&W=0/L=0,Age<=45,married) Pr(DI=1&W=0/L=1,Age<=45,married) Pr(DI=1&W=0/L=2,Age<=45,married) ///
									Pr(DI=1&W=0/L=0,Age>45,married) Pr(DI=1&W=0/L=1,Age>45,married) Pr(DI=1&W=0/L=2,Age>45,married) ///
									Pr(DI=1&W=1/L=0,Age<=45,single) Pr(DI=1&W=1/L=1,Age<=45,single) Pr(DI=1&W=1/L=2,Age<=45,single) ///
									Pr(DI=1&W=1/L=0,Age>45,single) Pr(DI=1&W=1/L=1,Age>45,single) Pr(DI=1/L=2&W=1,Age>45,single) ///
									Pr(DI=1&W=1/L=0,Age<=45,married) Pr(DI=1&W=1/L=1,Age<=45,married) Pr(DI=1&W=1/L=2,Age<=45,married) ///
									Pr(DI=1&W=1/L=0,Age>45,married) Pr(DI=1&W=1/L=1,Age>45,married) Pr(DI=1&W=1/L=2,Age>45,married) 
matrix colnames stockDI_def1 = EST SE	

***Table 4 Panel C
matlist stockDI_def1,cspec(o4& %25s | %5.4f & %5.4f & ) rspec(&-&&&&&&&&&&&&&&&&&&&&&&&&) title(DI rate(Year t),DS(Intvw. time t+1))
mat2txt, matrix(stockDI_def1) saving("$out/stockDImoments_work`1'") replace

********Stocks on DI by health, productive type, and age

matrix b=0
matrix V=0
cap program drop doit
program def doit
		 local ty = 1
		 levelsof type_fix, local(types)
		 foreach `ty' of local types {
         local i=0
	   while `i' < 2 {
		   local j=0	
		   while `j' < 3 {
				reg DI if DS==`j' & old==`i' & type_fix==`ty', vce(cluster newid)
				matrix b=b\e(b)
				matrix V=V\e(V)
	         local j=`j'+1
                 }
         local i=`i'+1
                 }
				 }
end
qui doit
matrix VV=0
forvalues j=2(1)19	{
	matrix temp=sqrt(V[`j',1])
	matrix VV=VV\temp
	}
matrix stockDI_def1=b[2..19,1],VV[2..19,1]


matrix rownames stockDI_def1 =	Pr(DI=1/L=0,Age<=45,T=1) Pr(DI=1/L=1,Age<=45,T=1) Pr(DI=1/L=2,Age<=45,T=1) ///
									Pr(DI=1/L=0,Age>45,T=1) Pr(DI=1/L=1,Age>45,T=1) Pr(DI=1/L=2,Age>45,T=1) ///
									Pr(DI=1/L=0,Age<=45,T=2) Pr(DI=1/L=1,Age<=45,T=2) Pr(DI=1/L=2,Age<=45,T=2) ///
									Pr(DI=1/L=0,Age>45,T=2) Pr(DI=1/L=1,Age>45,T=2) Pr(DI=1/L=2,Age>45,T=2) ///
									Pr(DI=1/L=0,Age<=45,T=3) Pr(DI=1/L=1,Age<=45,T=3) Pr(DI=1/L=2,Age<=45,T=3) ///
									Pr(DI=1/L=0,Age>45,T=3) Pr(DI=1/L=1,Age>45,T=3) Pr(DI=1/L=2,Age>45,T=3) 
									
matrix colnames stockDI_def1 = EST SE	

***Table 4 Panel C
matlist stockDI_def1,cspec(o4& %25s | %5.4f & %5.4f & ) rspec(&-&&&&&&&&&&&&&&&&&&) title(DI rate(Year t),DS(Intvw. time t+1))
mat2txt, matrix(stockDI_def1) saving("$out/stockDImoments_type`1'") replace


****Stocks on DI, by marital status, health, type, and age
********Stocks on DI

matrix b=0
matrix V=0
cap program drop doit
program def doit
local ty = 1
levelsof type_fix, local(types)
		 foreach `ty' of local types {
		 local m = 0
		 while `m' < 2 {
         local i=0
	   while `i' < 2 {
		   local j=0	
		   while `j' < 3 {
				reg DI if DS==`j' & old==`i' & married==`m' & type_fix == `ty', vce(cluster newid)
				matrix b=b\e(b)
				matrix V=V\e(V)
	         local j=`j'+1
                 }
         local i=`i'+1
                 }
				 local m = `m' + 1
				 }
				 }
end
qui doit
matrix VV=0
forvalues j=2(1)37	{
	matrix temp=sqrt(V[`j',1])
	matrix VV=VV\temp
	}
matrix stockDI_def1=b[2..37,1],VV[2..37,1]


matrix rownames stockDI_def1 =	Pr(DI=1/L=0,Age<=45,T=1,single) Pr(DI=1/L=1,Age<=45,T=1,single) Pr(DI=1/L=2,Age<=45,T=1,single) ///
									Pr(DI=1/L=0,Age>45,T=1,single) Pr(DI=1/L=1,Age>45,T=1,single) Pr(DI=1/L=2,Age>45,T=1,single) ///
									Pr(DI=1/L=0,Age<=45,T=1,married) Pr(DI=1/L=1,Age<=45,T=1,married) Pr(DI=1/L=2,Age<=45,T=1,married) ///
									Pr(DI=1/L=0,Age>45,T=1,married) Pr(DI=1/L=1,Age>45,T=1,married) Pr(DI=1/L=2,Age>45,T=1,married) ///
									Pr(DI=1/L=0,Age<=45,T=2,single) Pr(DI=1/L=1,Age<=45,T=2,single) Pr(DI=1/L=2,Age<=45,T=2,single) ///
									Pr(DI=1/L=0,Age>45,T=2,single) Pr(DI=1/L=1,Age>45,T=2,single) Pr(DI=1/L=2,Age>45,T=2,single) ///
									Pr(DI=1/L=0,Age<=45,T=2,married) Pr(DI=1/L=1,Age<=45,T=2,married) Pr(DI=1/L=2,Age<=45,T=2,married) ///
									Pr(DI=1/L=0,Age>45,T=2,married) Pr(DI=1/L=1,Age>45,T=2,married) Pr(DI=1/L=2,Age>45,T=2,married) ///
									Pr(DI=1/L=0,Age<=45,T=3,single) Pr(DI=1/L=1,Age<=45,T=3,single) Pr(DI=1/L=2,Age<=45,T=3,single) ///
									Pr(DI=1/L=0,Age>45,T=3,single) Pr(DI=1/L=1,Age>45,T=3,single) Pr(DI=1/L=2,Age>45,T=3,single) ///
									Pr(DI=1/L=0,Age<=45,T=3,married) Pr(DI=1/L=1,Age<=45,T=3,married) Pr(DI=1/L=2,Age<=45,T=3,married) ///
									Pr(DI=1/L=0,Age>45,T=3,married) Pr(DI=1/L=1,Age>45,T=3,married) Pr(DI=1/L=2,Age>45,T=3,married) 
matrix colnames stockDI_def1 = EST SE	

***Table 4 Panel C
matlist stockDI_def1, cspec(o4& %32s | %5.4f & %5.4f & ) rspec(&-&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&) title(DI rate(Year t),DS(Intvw. time t+1)) 
mat2txt, matrix(stockDI_def1) saving("$out/stockDImoments_martype`1'") replace



tab DS, gen(DS)

*keep if insample_head == 1
*keep if everwork==1
*restricting to everwork fixes health differences!
*but awkward... imposes different restrictions on old vs young
*use another measure of health (child health)
*within-type vs. across-type differences in health over marriage:
mat def marDS_across = J(6,2,0)
forvalues o=0(1)1{
local N1 0
local N0 0
forvalues j = 0(1)1{
forvalues i=1(1)3{
*total count of (mar, age) group size
quietly count if married==`j' & type_fix==`i' & old == `o'
local N`j' = `N`j''+`r(N)'
}

forvalues h=1(1)3{
forvalues i=1(1)3{
*number of (mar, age)  in type m
quietly count if married==`j' & type_fix==`i' & old == `o'
local num`h' =`r(N)'
quietly sum DS`h' if type_fix==`i' & old == `o'
mat marDS_across[`h'+3*`o',`j'+1] = marDS_across[`h'+3*`o',`j'+1] + (`num`h''*`r(mean)')/`N`j''
}
}
}
}

mat def marDS_total = J(6,2,0)
forvalues j = 0(1)1{
forvalues o=0(1)1{
forvalues h=1(1)3{
quietly sum DS`h' if married==`j' & type_fix != . & old==`o'
mat marDS_total[`h'+3*`o',`j'+1]=`r(mean)'
}
}
}
mat list marDS_across
mat list marDS_total

levelsof type_fix, local(types)
local ntypes: word count `types'
mat def marDS_all = J(6,`ntypes'*2+4,0)
local marcolnames 
foreach ty of local types{
local marcolnames `marcolnames' T`ty'_Single T`ty'_Married
}
local marcolnames `marcolnames' Weighted_Single Weighted_Married Raw_Single Raw_Married
*local marcolnames `""Low_Single" "Low_Married" "Med_Single" "Med_Married" "High_Single" "High_Married" "Weighted_Single" "Weighted_Married" "Raw_Single" "Raw_Married""'
local marrownames young_healthy young_mod young_sev old_healthy old_mod old_sev
foreach ty of local types{
forvalues j = 0(1)1{
forvalues o=0(1)1{
forvalues h=1(1)3{
quietly sum DS`h' if married==`j' & type_fix == `ty' & old==`o'
mat marDS_all[`h'+3*`o',(`j'+1)+2*(`ty'-1)]=`r(mean)'
}
}
}
}

forvalues j = 0(1)1{
forvalues o=0(1)1{
forvalues h=1(1)3{
mat marDS_all[`h'+3*`o',(`j'+1)+2*(3)]=marDS_across[`h'+3*`o',`j'+1]
mat marDS_all[`h'+3*`o',(`j'+1)+2*(4)]=marDS_total[`h'+3*`o',`j'+1]
}
}
}
mat colnames marDS_all= `marcolnames'
mat rownames marDS_all = `marrownames'

outtable using "$out/DIbymar`1'", mat(marDS_all)  replace asis

********Stocks on L
gen nodis		=DS		==0

matrix b=0
matrix V=0
reg sev if DI==1 & old==0, vce(cluster newid)
matrix b=b\e(b)
matrix V=V\e(V)
reg sev if DI==1 & old==1, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg mod if DI==1 & old==0, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg mod if DI==1 & old==1, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg nodis if DI==1 & old==0, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg nodis if DI==1 & old==1, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)

matrix VV=0
forvalues j=2(1)7	{
	matrix temp=sqrt(V[`j',1])
	matrix VV=VV\temp
	}
matrix stockL_def1=b[2..7,1],VV[2..7,1]


matrix rownames stockL_def1 =	Pr(L=2/DI=1,Age<=45) Pr(L=2/DI=1,Age>45)	///
									Pr(L=1/DI=1,Age<=45) Pr(L=1/DI=1,Age>45)	///
									Pr(L=0/DI=1,Age<=45) Pr(L=0/DI=1,Age>45)	///

matrix colnames stockL_def1 = EST SE	

***Table 4 Panel D
matlist stockL_def1,cspec(o4& %25s | %5.4f & %5.4f & ) rspec(&-&&&&&&) title(L propor.(Intvw. time t+1),DI(Year t))
mat2txt, matrix(stockL_def1) saving("$out/compDImoments`1'_agg") replace


matrix b=0
matrix V=0
reg sev if DI==1 & old==0 & married == 0, vce(cluster newid)
matrix b=b\e(b)
matrix V=V\e(V)
reg sev if DI==1 & old==1 & married == 0, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg mod if DI==1 & old==0 & married == 0, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg mod if DI==1 & old==1 & married == 0, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg nodis if DI==1 & old==0 & married == 0, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg nodis if DI==1 & old==1 & married == 0, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)

reg sev if DI==1 & old==0 & married == 1, vce(cluster newid)
matrix b=b\e(b)
matrix V=V\e(V)
reg sev if DI==1 & old==1 & married == 1, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg mod if DI==1 & old==0 & married == 1, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg mod if DI==1 & old==1 & married == 1, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg nodis if DI==1 & old==0 & married == 1, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg nodis if DI==1 & old==1 & married == 1, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)

matrix VV=0
forvalues j=2(1)13	{
	matrix temp=sqrt(V[`j',1])
	matrix VV=VV\temp
	}
matrix stockL_def1=b[2..13,1],VV[2..13,1]


matrix rownames stockL_def1 =	Pr(L=2/DI=1,Age<=45,single) Pr(L=2/DI=1,Age>45,single)	///
									Pr(L=1/DI=1,Age<=45,single) Pr(L=1/DI=1,Age>45,single)	///
									Pr(L=0/DI=1,Age<=45,single) Pr(L=0/DI=1,Age>45,single)	///
									Pr(L=2/DI=1,Age<=45,married) Pr(L=2/DI=1,Age>45,married)	///
									Pr(L=1/DI=1,Age<=45,married) Pr(L=1/DI=1,Age>45,married)	///
									Pr(L=0/DI=1,Age<=45,married) Pr(L=0/DI=1,Age>45,married)	

matrix colnames stockL_def1 = EST SE	

***Table 4 Panel D
matlist stockL_def1,cspec(o4& %25s | %5.4f & %5.4f & ) rspec(&-&&&&&&&&&&&&) title(L propor.(Intvw. time t+1),DI(Year t))
mat2txt, matrix(stockL_def1) saving("$out/compDImoments`1'") replace


********Stocks on L, by employment

matrix b=0
matrix V=0
gen sev0=sev==1&empnow==0
gen sev1=sev==1&empnow==1
gen mod0=mod==1&empnow==0
gen mod1=mod==1&empnow==1
gen nodis0=nodis==1&empnow==0
gen nodis1=nodis==1&empnow==1

forvalues l = 0(1)1{
reg sev`l' if DI==1 & old==0 & married == 0 , vce(cluster newid)
matrix b=b\e(b)
matrix V=V\e(V)
reg sev`l' if DI==1 & old==1 & married == 0, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg mod`l' if DI==1 & old==0 & married == 0, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg mod`l' if DI==1 & old==1 & married == 0, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg nodis`l' if DI==1 & old==0 & married == 0, vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg nodis`l' if DI==1 & old==1 & married == 0 , vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)

reg sev`l' if DI==1 & old==0 & married == 1 , vce(cluster newid)
matrix b=b\e(b)
matrix V=V\e(V)
reg sev`l' if DI==1 & old==1 & married == 1 , vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg mod`l' if DI==1 & old==0 & married == 1 , vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg mod`l' if DI==1 & old==1 & married == 1 , vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg nodis`l' if DI==1 & old==0 & married == 1 , vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
reg nodis`l' if DI==1 & old==1 & married == 1 , vce(cluster newid) 
matrix b=b\e(b)
matrix V=V\e(V)
}


matrix VV=0
forvalues j=2(1)25	{
	matrix temp=sqrt(V[`j',1])
	matrix VV=VV\temp
	}
matrix stockL_def1=b[2..25,1],VV[2..25,1]


matrix rownames stockL_def1 =	Pr(L=2&W=0/DI=1,Age<=45,single) Pr(L=2&W=0/DI=1,Age>45,single)	///
									Pr(L=1&W=0/DI=1,Age<=45,single) Pr(L=1&W=0/DI=1,Age>45,single)	///
									Pr(L=0&W=0/DI=1,Age<=45,single) Pr(L=0&W=0/DI=1,Age>45,single)	///
									Pr(L=2&W=0/DI=1,Age<=45,married) Pr(L=2&W=0/DI=1,Age>45,married)	///
									Pr(L=1&W=0/DI=1,Age<=45,married) Pr(L=1&W=0/DI=1,Age>45,married)	///
									Pr(L=0&W=0/DI=1,Age<=45,married) Pr(L=0&W=0/DI=1,Age>45,married)	///
									Pr(L=2&W=1/DI=1,Age<=45,single) Pr(L=2&W=1/DI=1,Age>45,single)	///
									Pr(L=1&W=1/DI=1,Age<=45,single) Pr(L=1&W=1/DI=1,Age>45,single)	///
									Pr(L=0&W=1/DI=1,Age<=45,single) Pr(L=0&W=1/DI=1,Age>45,single)	///
									Pr(L=2&W=1/DI=1,Age<=45,married) Pr(L=2&W=1/DI=1,Age>45,married)	///
									Pr(L=1&W=1/DI=1,Age<=45,married) Pr(L=1&W=1/DI=1,Age>45,married)	///
									Pr(L=0&W=1/DI=1,Age<=45,married) Pr(L=0&W=1/DI=1,Age>45,married)

matrix colnames stockL_def1 = EST SE	

***Table 4 Panel D
matlist stockL_def1,cspec(o4& %25s | %5.4f & %5.4f & ) rspec(&-&&&&&&&&&&&&&&&&&&&&&&&&) title(L propor.(Intvw. time t+1),DI(Year t))
mat2txt, matrix(stockL_def1) saving("$out/compDImoments_work`1'") replace




**************Flows onto DI
xtset newid year
g lag1DI=l1.DI
g lag2DI=l2.DI
g lag1DS=l1.DS
g lag2DS=l2.DS

matrix b=0
matrix V=0
cap program drop doit
program def doit
         local i=0
	   while `i' < 2 {
		   local j=0	
		   while `j' < 3 {
				reg DI if DS==`j' & ((lag2DI==0 & lag1DI == .) | lag1DI == 0) & old==`i'
				matrix b=b\e(b)
				matrix V=V\e(V)
	         local j=`j'+1
                 }
         local i=`i'+1
                 }
end
qui doit
matrix VV=0
forvalues j=2(1)7	{
	matrix temp=sqrt(V[`j',1])
	matrix VV=VV\temp
	}
matrix flowsinto_2year=(b[2..4,1]\b[5..7,1]),(VV[2..4,1]\VV[5..7,1])


matrix rownames flowsinto_2year =			Pr(DI=1/LagDI=0,L=0,Age<=45) Pr(DI=1/LagDI=0,L=1,Age<=45) Pr(DI=1/LagDI=0,L=2,Age<=45) ///
											Pr(DI=1/LagDI=0,L=0,Age>45) Pr(DI=1/LagDI=0,L=1,Age>45) Pr(DI=1/LagDI=0,L=2,Age>45) 
matrix colnames flowsinto_2year = EST SE	

***Table 4 Panel E
matlist flowsinto_2year,cspec(o4& %28s | %5.4f & %5.4f & ) rspec(&-&&&&&&) title(Pr(DI(t)=1|DI(t-2)=0,L(t+1),Age))
mat2txt, matrix(flowsinto_2year) saving("$out/flowDImoments`1'_agg") replace


***************FLOWS ONTO DI BY MARITAL STATUS
xtset newid year

matrix b=0
matrix V=0
cap program drop doit
program def doit
		 local m = 0
		 while `m' < 2 {
         local i=0
	   while `i' < 2 {
		   local j=0	
		   while `j' < 3 {
				reg DI if DS==`j' & ((lag2DI==0 & lag1DI == .) | lag1DI == 0) & old==`i' & married == `m' 
				matrix b=b\e(b)
				matrix V=V\e(V)
	         local j=`j'+1
                 }
         local i=`i'+1
                 }
				 local m = `m'+1
				 }
end
qui doit
matrix VV=0
forvalues j=2(1)13	{
	matrix temp=sqrt(V[`j',1])
	matrix VV=VV\temp
	}
matrix flowsinto_2year=(b[2..4,1]\b[5..7,1]\b[8..10,1]\b[11..13,1]),(VV[2..4,1]\VV[5..7,1]\VV[8..10,1]\VV[11..13,1])


matrix rownames flowsinto_2year =			Pr(DI=1/LagDI=0,L=0,Age<=45,M=0) Pr(DI=1/LagDI=0,L=1,Age<=45,M=0) Pr(DI=1/LagDI=0,L=2,Age<=45,M=0) ///
											Pr(DI=1/LagDI=0,L=0,Age>45,M=0) Pr(DI=1/LagDI=0,L=1,Age>45,M=0) Pr(DI=1/LagDI=0,L=2,Age>45,M=0) ///
											Pr(DI=1/LagDI=0,L=0,Age<=45,M=1) Pr(DI=1/LagDI=0,L=1,Age<=45,M=1) Pr(DI=1/LagDI=0,L=2,Age<=45,M=1) ///
											Pr(DI=1/LagDI=0,L=0,Age>45,M=1) Pr(DI=1/LagDI=0,L=1,Age>45,M=1) Pr(DI=1/LagDI=0,L=2,Age>45,M=1)
matrix colnames flowsinto_2year = EST SE	

***Table 4 Panel E
matlist flowsinto_2year,cspec(o4& %28s | %5.4f & %5.4f & ) rspec(&-&&&&&&&&&&&&) title(Pr(DI(t)=1|DI(t-2)=0,L(t+1),Age))
mat2txt, matrix(flowsinto_2year) saving("$out/flowDImomentsbymar`1'") replace

*************** FLOWS FOR CONTINUING DI BENEFICIARIES 
xtset newid year

matrix b=0
matrix V=0
cap program drop doit
program def doit
         local i=0
	   while `i' < 2 {
		   local j=0	
		   while `j' < 3 {
				reg DI if  DS==`j' & ((lag2DI==1 & lag1DI == .) | lag1DI == 1) & old==`i' ,vce(cluster newid) 
				matrix b=b\e(b)
				matrix V=V\e(V)
	         local j=`j'+1
                 }
         local i=`i'+1
                 }
end
qui doit
matrix VV=0
forvalues j=2(1)7	{
	matrix temp=sqrt(V[`j',1])
	matrix VV=VV\temp
	}
matrix flowsinto_2year=(b[2..4,1]\b[5..7,1]),(VV[2..4,1]\VV[5..7,1])


matrix rownames flowsinto_2year =			Pr(DI=1/LagDI=1,L=0,Age<=45) Pr(DI=1/LagDI=1,L=1,Age<=45) Pr(DI=1/LagDI=1,L=2,Age<=45) ///
											Pr(DI=1/LagDI=1,L=0,Age>45) Pr(DI=1/LagDI=1,L=1,Age>45) Pr(DI=1/LagDI=1,L=2,Age>45) 
matrix colnames flowsinto_2year = EST SE	

***Table 4 Panel E
matlist flowsinto_2year,cspec(o4& %28s | %5.4f & %5.4f & ) rspec(&-&&&&&&) title(Pr(DI(t)=1|DI(t-2)=0,L(t+1),Age))
mat2txt, matrix(flowsinto_2year) saving("$out/contflowDImoments`1'_agg") replace



*************** FLOWS FOR CONTINUING DI BENEFICIARIES BY MARITAL STATUS
xtset newid year

matrix b=0
matrix V=0
cap program drop doit
program def doit
		 local m = 0
		 while `m' < 2 {
         local i=0
	   while `i' < 2 {
		   local j=0	
		   while `j' < 3 {
				reg DI if DS==`j' & ((lag2DI==1 & lag1DI == .) | lag1DI == 1) & old==`i' & married == `m' ,vce(cluster newid) 
				matrix b=b\e(b)
				matrix V=V\e(V)
	         local j=`j'+1
                 }
         local i=`i'+1
                 }
				 local m = `m'+1
				 }
end
qui doit
matrix VV=0
forvalues j=2(1)13	{
	matrix temp=sqrt(V[`j',1])
	matrix VV=VV\temp
	}
matrix flowsinto_2year=(b[2..4,1]\b[5..7,1]\b[8..10,1]\b[11..13,1]),(VV[2..4,1]\VV[5..7,1]\VV[8..10,1]\VV[11..13,1])


matrix rownames flowsinto_2year =			Pr(DI=1/LagDI=1,L=0,Age<=45,M=0) Pr(DI=1/LagDI=1,L=1,Age<=45,M=0) Pr(DI=1/LagDI=1,L=2,Age<=45,M=0) ///
											Pr(DI=1/LagDI=1,L=0,Age>45,M=0) Pr(DI=1/LagDI=1,L=1,Age>45,M=0) Pr(DI=1/LagDI=1,L=2,Age>45,M=0) ///
											Pr(DI=1/LagDI=1,L=0,Age<=45,M=1) Pr(DI=1/LagDI=1,L=1,Age<=45,M=1) Pr(DI=1/LagDI=1,L=2,Age<=45,M=1) ///
											Pr(DI=1/LagDI=1,L=0,Age>45,M=1) Pr(DI=1/LagDI=1,L=1,Age>45,M=1) Pr(DI=1/LagDI=1,L=2,Age>45,M=1)
matrix colnames flowsinto_2year = EST SE	

***Table 4 Panel E
matlist flowsinto_2year,cspec(o4& %28s | %5.4f & %5.4f & ) rspec(&-&&&&&&&&&&&&) title(Pr(DI(t)=1|DI(t-2)=0,L(t+1),Age))
mat2txt, matrix(flowsinto_2year) saving("$out/contflowDImomentsbymar`1'") replace



**************Flows off DI
gen noDI=1-DI

forvalues i=0(1)2{
matrix b=0
matrix V=0
tab married if ((lag2DI==1 & lag1DI == .) | lag1DI == 1) & old==0 & DS == `i'
reg noDI married if ((lag2DI==1 & lag1DI == .) | lag1DI == 1) & old==0 & DS == `i',vce(cluster newid) 
tab married if ((lag2DI==1 & lag1DI == .) | lag1DI == 1) & old==1 & DS == `i'
reg noDI married if ((lag2DI==1 & lag1DI == .) | lag1DI == 1) & old==1 & DS == `i',vce(cluster newid) 
}
*If anything, married people leave DI MORE OFTEN.

/*
matrix VV=0
forvalues j=2(1)3	{
	matrix temp=sqrt(V[`j',1])
	matrix VV=VV\temp
	}
matrix flowsoff=(b[2..3,1]),(VV[2..3,1])
matrix rownames flowsoff =	Pr(DI=0/DI'=1,Age<=45) Pr(DI=0/DI'=1,Age>45) 
matrix colnames flowsoff = EST SE	
***Table 8
matlist flowsoff,cspec(o4& %28s | %5.4f & %5.4f & ) rspec(&-&-) title(Pr(DI(t)=0|DI(t-j)=1,Age))
*/
************Asset distribution by age
preserve
forvalues p=1(1)99{
local coll (p`p') asset_p`p'=w_with_h `coll'
}
collapse `coll' , by(age type_fix)

forvalues p=1(1)99{
gen asset_p`p'_smooth =.
foreach typ of local types{
lowess asset_p`p' age if type_fix==`typ', gen(asset_p`p'_smooth_`typ')
replace asset_p`p'_smooth = asset_p`p'_smooth_`typ' if type_fix==`typ'
drop asset_p`p'_smooth_`typ'
}
}
keep age type_fix asset*smooth
foreach x of varlist asset*smooth{
replace `x' = 0 if `x' < 0
}
rename asset*_smooth asset_smooth*
reshape long asset_smooth_p, i(age type_fix) j(p)
rename asset_smooth_p asset

replace asset = 0 if asset < 0
export delimited using "$out/assetdist_byage`1'.csv", replace

macro shift
restore

clear
 
